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ABSTRACT 


Data from one-dimensional (spherically symmetric) analyses was used to 
examine the effects of compressibility and gas energy on the dynamic behavior of an 
explosion gas bubble, by comparing the bubble's behavior with experimental results 
and with analytical results which neglect these factors. Results from two-dimensional 
(axially symmetric) analyses were used to investigate the behavior of a deep 
explosion gas bubble in the vicinity of plane rigid or constant pressure boundaries. 
Previous analytical research into explosion gas bubbles near such boundaries has 
primarily led to results of a qualitative nature, owing to a complete breakdown of the 
assumptions made in the analysis at the critical juncture. In the present investigation, 
it was found possible to characterize the effect of the boundary surface on both the 
change in the first oscillation period of the bubble and its location at the end of the 
first oscillation cycle. For a broad range of bubble - boundary standoff distances, 
these semi-empirical characterizations have a functional form particularly suitable for 
extension of the quantitative results of this investigation to other explosive charge 
types, weights, and depths, as has been done for the Willis formula for the free-field 


oscillation period of explosion gas bubbles. 
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I, INTRODUCTION 


The survivability of naval ships and submarines is of supreme importance to those who 
design, build and sail them. Underwater explosions, created by the detonation of mines or 
torpedoes near naval vessels, clearly represent a significant threat to that survivability. The 
capacity of such weapons to cause major damage to naval vessels has been recognized since 
the beginning of the age of modern naval warfare, and has been the subject of extensive 
investigation since the earliest years of the Second World War. 

One result of the allied experience during the Second World War was the realization 
that the pressure wave resulting from the pulsation (oscillation) of the bubble of gases 
produced by the detonation of an explosive material underwater can be as important as the 
primary shock wave in causing damage to a naval vessel. As Hicks (1970a) has noted, this 
bubble pulsation is the primary cause of the observed "whipping" of ships, and often has a 
period near that of the hull girder's fundamental flexural vibration mode. 

This will, in the worst case, break the back of the ship. Even in less severe situations, 
the oscillation of the ship's hull may be significantly increased by the timing of the arrival of 
the secondary pressure pulse radiated when the bubble collapses to its first minimum volume. 
This can cause appreciable damage to sensitive internal equipment. Furthermore, the 
increased use of commercial off-the-shelf (COTS) equipment within modern ships, while 
fiscally necessary and undoubtedly beneficial in the long run, may increase the damage 
potential of explosion gas bubble oscillations. 


Fundamental analysis shows that low frequency vibration input can significantly 


increase the resulting response in certain underdamped system configurations. Equipment not 
designed to withstand the increased response could easily suffer damage that, depending upon 
ship design, reduces the war fighting capability of the ship. In an age of over-the-horizon 
saturation targeting, such a calamity may be as bad as breaking the back of the ship. 

A significant investment of manpower and time has been expended in analyzing the 
phenomena associated with the pulsation of explosion gas bubbles, both during and since the 
Second World War. However, these studies have, for the most part, relied upon many 
simplifying assumptions. The validity of some of these assumptions were known beforehand 
to be very inaccurate during certain phases of an explosion gas bubble's motion. In fact, some 
of these assumptions are least accurate during the times when key explosion gas bubble 
phenomena are occurring. Predictions of explosion gas bubble behavior from these studies 
are therefore more likely to be qualitatively correct than quantitatively accurate. This is an 
important consideration because ship design decisions are based, in part, upon these 
predictions. 

The next section describes some of the important phenomena associated with 
underwater explosion gas bubbles. Following this is a synopsis of the more significant 
research that has been previously conducted into underwater explosion gas bubbles. The third 
and final section of this introductory chapter discusses the objectives of this present research 
effort, and presents an overview of the remaining chapters. 

A. BUBBLE PHENOMENA 
Consider a spherical explosive charge located underwater at a considerable distance 


from any boundary surface. When this charge is detonated, the explosive material is very 


rapidly converted into high pressure gaseous reaction products, and a high pressure shock 
wave is transmitted to and propagated through the surrounding fluid. About 50% of the 
initial chemical energy of the explosive is transmitted to the fluid in this initial shock wave 
(Snay, 1957). 

Because this shock wave propagates radially outward, the amplitude of the shock 
wave decreases with increasing distance from the center of the charge. The fluid behind the 
shock wave front attains a large outward radial velocity. This causes the pressure at some 
distance behind the shock wave front to drop below the hydrostatic pressure of the fluid 
ahead of this wave front. 

The pressure in the gas bubble is significantly reduced after emission of the primary 
shock wave, but is still significantly higher than the hydrostatic pressure in the surrounding 
fluid. This, in conjunction with the radial outflow of fluid from the vicinity of the bubble, 
causes the bubble to expand rapidly. This expansion continues for a relatively long time, the 
pressure in the bubble decreasing as the bubble volume increases. The pressure in the gas 
bubble eventually falls below hydrostatic pressure, but the expansion persists because of the 
inertia of the outward flowing fluid (Cole, 1948, p. 8). 

Eventually, the fluid immediately surrounding the bubble comes to rest. This, together 
with the very low pressure in the bubble at this time, causes the bubble to begin contracting. 
The fluid in the immediate vicinity of the bubble begins flowing inward, accelerating the 
contraction of the bubble. As the volume of the bubble decreases, its internal pressure 
increases. The pressure within the gas bubble eventually becomes significantly higher than 


the hydrostatic pressure of the surrounding fluid, and the contraction of the bubble is reversed 


abruptly. The inertia of the surrounding fluid together with the compressibility of the bubble 
gases and the surrounding fluid "thus provide the necessary conditions for an oscillating 
system, and the bubble does in fact undergo repeated cycles of expansion and contraction" 
(Cole, 1948, p. 8). 

Figure 1-1 illustrates the bubble radius vs. time behavior typically observed in an 
underwater explosion. Note that the maximum radius of the bubble decreases during 
subsequent expansions. This is due, at least in part, to the emission of secondary pressure 
pulses when the bubble radius is near its minimum values. About 66% of the energy 
remaining in the gas bubble after emission of the primary shock wave is lost during the first 
expansion-contraction cycle (Cole, 1948, p. 283). Successive bubble pulses are thus weaker, 
and generally only the first pulse is of practical significance (Cole, 1948, p. 10). 

Figure 1-2 shows a typical pressure vs. time curve at fixed distance from an 
underwater explosion. Although the peak pressure in the secondary pressure pulse is much 
less than the peak shock wave pressure, the duration of this pulse 1s greater. The impulses 
due to the primary shock wave and the first of the bubble pulses are typically comparable 
(Cole, 1948, p. 364). 

Additional phenomena are observed when an explosion occurs near a boundary 
surface, such as the surface or bottom of the ocean or a nearby structure. When an explosion 
occurs near a relatively ngid boundary such as the hull of a naval vessel or a hard ocean 
bottom, the explosion gas bubble migrates toward the boundary during its contraction. The 
opposite effect is seen when an explosion occurs near a free surface (Snay, 1957). The shape 


of migrating gas bubbles has also been seen to be non-spherical, as illustrated in Figure 1-3. 
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Figure 1-1. Typical Bubble Radius vs. Time Curve [From Swift and Decius (1950)] 
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Figure 1-3. Typical Shape for a Migrating Underwater Explosion Gas Bubble [From 
Cole (1948)] 


Since the secondary pressure pulse is emitted when a bubble is near its minimum 
volume, the migration of a bubble determines where this bubble pulse is emitted from. For 
an explosion occurring near a naval ship, the migration of the bubble can cause the bubble 
pulse to occur much closer to the hull, potentially increasing the damage caused to the ship 
(Hicks, 1970b). 

The presence of boundary surfaces also affects the period of the oscillation of the 
bubble. A rigid boundary causes the period of oscillation of the bubble to increase, and a free 
surface causes it to decrease (Herring, 1950). Thus, if an explosion occurs near the hull of 
a ship, the presence of the ship itself can cause the time interval between arrival of the primary 
shock wave and secondary pressure pulse to become either closer to or further from the 
fundamental period of the hull girder, significantly affecting the amount of whipping 
undergone by the ship. 

B. PREVIOUS RESEARCH INTO EXPLOSION GAS BUBBLES 

Perhaps the earliest contribution related directly to an appreciation and understanding 
of bubbles was made by Reynolds (1894), who noted the formation of vapor cavitation 
bubbles in water flowing through a constricted pipe. The first analysis of the dynamic 
behavior of bubbles was made by Rayleigh (1917), who used conservation of momentum to 
derive an equation for the collapse of a spherical void. 

Lamb (1923) studied the expansion phase of a gas bubble, and derived an expression 
relating the maximum radius of a bubble to the depth and total energy, neglecting fluid 
compressibility. Ramsauer (1923) conducted small scale experiments with guncotton which 


showed reasonably good agreement with Lamb's predicted relationship between maximum 


bubble radius and depth. The work of Lamb is the basis for the standard semi-empirical 
scaling formula for the maximum radius of an explosion gas bubble when buoyancy can be 
neglected and no boundary surfaces are near the charge. This situation is referred to as the 
"free-field" case. This scaling formula defines a relationship between the maximum radius of 
the bubble, the charge weight, and the charge depth. It is termed "semi-empirical" because 
it depends upon a proportionality constant which varies between different types of explosives, 
and must be determined through experimentation. 

Willis (1941) was the first to succeed in integrating the energy equation for the 
noncompressive radial motion of free-field explosion gas bubbles, and thereby derive an 
expression for the period of the first oscillation cycle of the bubble. His work is the basis for 
the standard semi-empirical formula for the free-field first bubble oscillation period of an 
explosion gas bubble. This formula defines a scaling relationship for different charge weights 
and depths, and depends upon an experimentally determined proportionality constant. Willis 
also conducted experiments to confirm the predicted relationship between charge size, depth, 
and bubble period, and found good agreement. This result has since been confirmed by other 
researchers, particularly as to the functional dependency of the first bubble period on the 
charge depth (Cole, 1948, pp. 280-281). 

Herring (1941, 1950) made a number of important contributions to the current 
understanding of explosion gas bubble phenomena. He predicted, based upon an approximate 
theoretical development, that one of the effects of fluid compressibility would be to make the 
bubble radius versus time curve for a free-field bubble asymmetrical, with the contraction 


being slower the expansion. He conducted an analysis indicating that the radiation of energy, 


if fluid compressibility was allowed for, would be negligible when the bubble is large, but that 
an appreciable radiation of energy could take place at the first bubble minimum. 

Herring's most important contribution was the advancement of an approximate theory 
to account for the effect of simple rigid boundaries or free surfaces (constant pressure 
boundaries) on explosion gas bubbles during the first oscillation cycle. In this theory, the 
effects of the boundary surface are treated as small perturbations of the motion of the gas 
bubble in the absence of the boundary. The boundary is assumed to be remote enough so that 
the standoff distance to the boundary is large compared to the maximum radius of the bubble, 
which is assumed to remain nearly spherical. This theory is thus only potentially valid when 
the boundary surface is fairly remote from the gas bubble. Furthermore, the theory neglects 
fluid compressibility entirely. 

One of the effects of fluid compressibility is a finite wave speed in the fluid; in an 
incompressible fluid, the wave velocity is infinite. Thus in a compressible fluid there 1s a finite 
distance beyond which the presence of a boundary cannot affect an explosion gas bubble 
during the first oscillation cycle, because even the primary shock wave resulting from 
detonation of the charge will not have reached the boundary by the end of the first cycle. The 
theory developed by Herring to account for simple boundary surfaces therefore cannot be 
correct at large standoff distances, as it assumes an incompressible fluid. 

Although Herring's theory for predicting the effect of these simple types of boundary 
surfaces is still in use, e.g. in the computer code MSWHIP (Hicks, 1971) (which 1s still being 
used by the United States Navy to predict the response of naval vessels to underwater 


explosion gas bubbles), the above discussion begs an answer to the question of when It is 
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valid. It was never intended to be accurate at small standoff distances, and because it neglects 
fluid compressibility, it is not valid at large standoff distances. 

Another important contnbution to the current understanding of explosion gas bubble 
phenomena was made by Taylor (1943), who denved an expression for the vertical migration 
of a spherically symmetric bubble due to gravity (buoyancy) in the absence of nearby 
boundary surfaces, neglecting fluid compressibility. In addition to neglecting fluid 
compressibility, Taylor's theory is based upon the assumption that artificial forces, which do 
no work, act to keep the gas bubble spherical. Actual experimental results typically show 
significant bubble shape departure from spherical, as illustrated in Figure 1-3. In the present 
application of Taylor's theory for the effects of buoyancy on explosion gas bubbles, use is also 
made of Herring's theory for the effect of the free surface (if a free surface is not nearby, the 
bubble is generally deep enough so that buoyancy is not significant). 

It is perhaps then not surprising that predictions of explosion gas bubble behavior from 
the combined Herring-Taylor theory do not agree very well with experimental results. For 
example, an experimental measurement by Bryant (1950) near a free surface showed a peak 
bubble migration velocity of 180 ft/s (55 m/s), as compared to a calculated value of 890 ft/s 
(271 m/s). 

Numerous researchers have conducted underwater explosion experiments. Notable 
among these is Swift and Decius (1950), who conducted a series of experiments to determine 
relatively accurate values for the proportionality constants in the semi-empirical formulas for 
the free-field maximum bubble radius and first bubble period, for various types of explosives. 


Chertock (1952) developed an approximate method to calculate the bubble pulse 


1] 


induced whipping response of naval ships and submarines, for explosions occurring very 
remote from the hull. Hicks (1970a, 1970b, 1971) determined that Chertock's methodology 
might be applied at smaller standoff distances with simple modifications, and then further 
modified this procedure to account for bubble migration, using the Herring-Taylor theory. 
In an attempt to account for the known problems with this theory, his procedure introduces 
an artificially large drag coefficient, the value of which he has chosen based upon very limited 
experimental data. This is justified as being "the best [solution] available until fresh 
experimental results (particularly for the near field) or a more rational theory become 
available" (Hicks, 1970b). The computer code developed by Hicks, MSWHIP, thus far 
appears to generally give reasonably accurate results for distant charges; the validity of its 
predictions at short and intermediate standoff distances is still largely unknown. 

Some recent efforts have examined methods to study underwater explosion gas bubble 
phenomena under laboratory conditions. Schmidt et. al (1987) have studied the problem of 
conducting very small scale explosive testing (charge weights of about 0.2 g) while 
maintaining complete similarity, by using a centrifuge to obtain gravitational accelerations of 
up to 500 times the standard terrestrial gravitational acceleration. Chahine et. al. (1995) have 
developed a procedure for generating plasma bubbles using a very high energy spark and 
scaling these plasma bubbles to represent explosion gas bubbles. 

As can be seen from the brief summary of previous work given above, much of the 
current understanding and theoretical development regarding underwater explosions was a 
result of allied research efforts during the Second World War. Many of the more important 


papers from this era can be found in the three volume work Underwater Explosion Research: 


IZ 


A Compendium of British and American Reports (Office of Naval Research, 1950). A 
comprehensive summary of the research from this period is provided by Cole (1948). 
CG OBJECTIVES OF CURRENT RESEARCH 

The above synopsis of previous research into explosion gas bubble phenomena points 
out several areas in which the existing understanding of these phenomena is inadequate. The 
existing theoretical understanding of the dynamic behavior of explosion gas bubbles is based 
upon neglecting fluid compressibility entirely. Only limited, qualitative estimates of what the 
affects of fluid compressibility might be have been made. And in addition to neglecting fluid 
compressibility, the existing theory descnbing the dynamic behavior of explosion gas bubbles 
near simple rigid or constant pressure boundary surfaces does not cover the situation in which 
the bubble is fairly close to the boundary. Even for the case in which the bubble is not too 
close to the boundary surface, the dynamic behavior of the bubble predicted using the existing 
theory is known to be more qualitatively correct than quantitatively accurate. 

The research described in this study was undertaken to rectify some of these 
shortcomings in the current understanding of explosion gas bubble phenomena. It was 
believed that, by using a modern numerical analysis computer program based upon a finite 
control volume (Eulerian) method, it would be possible to quantitatively investigate explosion 
gas bubble phenomena which had heretofore been understood only qualitatively. 

The first goal of this study was to examine the feasibility of using a finite control 
volume numerical analysis technique to investigate underwater explosion phenomena, and if 
this method proved serviceable, to investigate the effects of fluid compressibility and internal 


gas energy on the dynamic behavior of explosion gas bubbles. The second objective was to 


13 


investigate the effects of plane ngid and constant pressure boundaries on important aspects 
of the dynamic behavior of explosion gas bubbles, including fluid compressibility and internal 
energy, and to quantitatively characterize these effects. 

Two specific aspects of the dynamic behavior of bubbles near boundaries were 
investigated. The first of these was the change in the period of the bubble caused by the 
boundary. The bubble period 1s important because the time interval between arrival of the 
primary shock wave and the first bubble pulse can be very near the fundamental bending mode 
period of a ship's hull girder. Thus, if the boundary causes this time interval to become even 
closer to this hull girder bending mode period, the resultant whipping can be increased due 
to resonance. As discussed earlier, the impulses due to the pnmary shock wave and the first 
bubble pulse are comparable in magnitude, while the impulses from subsequent bubble pulses 
are much smaller. The other aspect of a bubble's dynamic behavior that was investigated was 
the migration of the bubble at the end of the first expansion-contraction cycle, either towards 
or away from the boundary. Because pressures fall off with increasing distance, it is 
important to know where the bubble is located when the first bubble pulse is emitted. 

The numenical analysis approach used for this research is described in the next chapter. 
Chapter III descnbes the application of this analysis approach to an investigation of the affects 
of fluid compressibility and internal gas energy on a free-field explosion gas bubble, and 
compares results from this analysis approach with both expenmental and analytical results. 
The effects of ngid and constant pressure boundaries on explosion gas bubbles, including fluid 
compressibility and internal gas energy, are examined in Chapter IV. Chapter V then 


summarizes the results from this study. 
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This dissertation illustrates the application of a finite volume based numerical analysis 
technique for investigation of underwater explosion gas bubbles. This technique was used to 
investigate the effects of fluid compressibility and internal energy on a free-field explosion gas 
bubble. This analysis procedure was also used to quantitatively characterize the effects of 
rigid and constant pressure boundaries on explosion gas bubbles, when internal energy and 
fluid compressibility are not neglected. The quantitative characterizations derived have a 
form suitable for use in scaling these results for other explosive charge types, weights, and 


depths. 


> 





H. NUMERICAL SOLUTION TECHNIQUE 


Because underwater explosions involve the flow of liquids and gases, and a high 
pressure shock wave, numerical analysis of underwater explosion phenomena in the region 
near the charge using standard Lagrangian based finite element programs is not practical. In 
a Lagrangian based finite element program, the elements deform in response to the pressure 
on the faces of the element. This leads to elements in the region of the shock front being 
severely deformed (crushed), and the time step size per iteration becoming very small, so that 
the solution advances very slowly in time. 

The numerical analyses described in this study were therefore conducted using an 
Eulerian based finite volume program. The computer program MSC/DYTRAN (The 
MacNeal Schwendler Corporation, 1995) was used for the numerical analyses conducted 
during this investigation. The analyses used the multi-material Eulerian processor in this 
program, which is based upon the computer program MSC/PISCES (The MacNeal 
Schwendler Corporation, 1991). This processor provides the flexibility of general 
connectivity for Eulerian finite volumes, so the finite volumes are not restricted as to shape. 
This feature is useful when a large volume of fluid must be discretized, but the area of interest 
is relatively small. Other Eulerian based programs typically require a uniform discretization, 
which, for analysis of an underwater explosion gas bubble, would require either an 
unacceptably coarse mesh or an unacceptably large number of finite volumes. This program 
also provides a constitutive equation suitable for modeling of explosive charge detonation, 


and a detonation wave front algorithm for detonating an explosive material. 
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This processor uses the control volume method to solve the basic conservation 
equations in space. The integral form of these equations, for conservation of mass, linear 


momentum, and total energy, are 


ee 
ral | J tame? ways -| J eed mards)+ J ac’ Ce" 
2H ant eaten 


where 7 is the stress tensor. For the hydrodynamic material models assumed for the 
numerical analyses of this investigation, 7 has the value -p along the main diagonal, and zero 
everywhere else. Viscosity and thermal conductivity are neglected in these equations. 

These equations are solved for each finite volume. A one point approximation is used 
to solve these equations (the value at the geometric center of the control volume), in 
conjunction with interpolated velocities and pressures at the faces of the finite volumes. 

A first order "Donor-Acceptor" scheme is used for material transport. Transported 
quantities are subtracted from donor cells and added to the acceptor cells, based upon the 
donor cell values and the velocity at the common face. The face velocity for the face 


connecting finite volumes m and 77 Is 
= b 
Urace 3 im mee ) (2-4) 


and, during time step df, the volume transport 1s 
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dV=i, -dSdt (2-5) 


face 


The normal component of the face velocity determines which cell is a donor and which is an 
acceptor (1.e. in which direction flow occur). If, for example, cell m is the donor, the mass, 
momentum, and energy transport from element mm to 7 is 


dM=p dV=p a, -dSdt (2-6) 


m face 


d(Mi)=p id dV=p i (i,.,-dS)dt O=7) 
de= p_¢,,aV= Pi? mit ‘ace Pout (2-8) 


The solution in time is computed with this program using an explicit central finite 
difference method. With this method, a new time step size is calculated after each forward 
time step. This new time step is calculated such that the Courant criterion for a stable 


solution, 


dt=S—— (2-9) 


u+c 


is satisfied, where S is a factor of safety, Z is the smallest element dimension, c is the acoustic 
wave velocity within the element, and w is the partical velocity within the element. This is 
calculated for all elements, and the smallest value obtained is used as the new time step. 


A complete time step algorithm consists of several phases. The momentum change 
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due to the impulse (from the face pressures) is calculated and used to update the velocity in 
each finite volume. Similarly, the work done by the pressure on each finite volume face is 
calculated and used to update the total energy for each finite volume. Mass, momentum and 
total energy transport across all faces are computed, and used to update these values in all 
finite volumes. The density and velocity in each finite volume are then updated (e.g. density 
equals new mass divided by fixed volume). The specific internal energy in each finite volume 
is computed based upon the new specific total energy and the new specific kinetic energy. 
And the pressure in each finite volume is updated using the new density and specific internal 
energy in the state equation. 

To define a multimaterial Eulerian model with this computer program, several inputs 
are required. In addition to defining the geometry of the mesh, constitutive equation data 
must be input for the different materials to be used. Initial conditions must be assigned to all 
of the finite volumes in the model; for a multimaterial Eulerian problem these would include 
specification of what matenals are initially in what finite volumes, and what the initial density 
and specific internal energy is in the different regions within the model. Any boundary 
conditions which are not to be left as the default "no flow" boundary condition must be 
specified. And, if the problem will use an explosive material, one or more detonation points, 
detonation velocities, and detonation times must be specified. Finally, the ending time for the 
analysis must be specified, and the program must be told precisely which output quantities are 


desired at what frequency. 
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Il. BASIC EXPLOSION GAS BUBBLE BEHAVIOR 


The effect of fluid compressibility and gas internal energy on the behavior of a free- 
field explosion gas bubble 1s investigated in this chapter. As discussed in the introduction, the 
basic charactenstics of an explosion gas bubble are a rapid expansion to a much larger volume 
than the initial volume of the charge, a prolonged interval of time in this expanded state, and 
a rapid collapse back to a small volume. This oscillation cycle then repeats, but with a 
decreasing period between collapses and a decrease in the maximum volume to which the 
bubble expands. 

In the first section in this chapter, the basic theoretical equations for the 
noncompressive motion of a free-field explosion gas bubble are developed. The second 
section describes the details of the numerical analysis models used in this investigation. 
Results from these analyses are presented in the last section. 

A. THEORETICAL DEVELOPMENT 

The classical expression for the motion of an explosion gas bubble in the absence of 
a nearby boundary surface neglects buoyancy resulting from pressure variation around the 
bubble and compressibility of the surrounding fluid. With these restrictions, the flow of fluid 


outside the bubble is radial, and conservation of energy gives 


3{ 4x 3){ dR > An 3 
—| — —| +—P.R°+E(R)=Y 3-] 
o{ $25,2°)| SB)’ «AE p,R? 20 G1) 


where R is the radius of the bubble, p, is the density in the surrounding fluid, Py is the 
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hydrostatic pressure in the surrounding fluid, E(R) is the internal energy of the gas bubble 
when it has radius R, and Y ts the total energy (a constant) (Cole, 1948, p.273). 

In this expression, the first term represents the kinetic energy of the flow in the 
surrounding fluid, and the second term is the work done against hydrostatic pressure in 
expanding the bubble to radius R. By neglecting the internal energy of the gas, which has a 
relatively small value over much of the oscillation cycle of an explosion gas bubble, the total 


energy can be expressed as 
4 
Y= P Rina (3-2) 


where &,,,. is the maximum radius of the bubble (Cole, 1948, p.274-275). This expression 
simply reflects the fact that the gas bubble must have a maximum radius when dR/dt is zero. 
Using equation (3-2) in equation (3-1), with E(R) taken as zero, equation (3-1) can be 


separated and integrated to give 


t=(3p,/2P,)!? ip *U(R,,,./a°-1]!?da (3-3) 
0 


where R, is the initial radius of the gas bubble. By taking R, as zero (the initial radius of the 
bubble is much smaller than the maximum radius) and transforming and integrating equation 


(3-3), Willis (1941) arrived at the result 


p 1/2 
ji ol 7 (3-4) 
Py 
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where T is the first oscillation period of the bubble. Substituting equation (3-2) into (3-4), 
and using the fact that the total energy Y is proportional to the weight W of the explosive 
charge gives 


wis 
P 5/6 


ane 





(3-5) 


where K; is a constant particular to a given type of explosive. Equation (3-2) can be re- 


arranged to give 


wis 
ae eT (3-6) 
0 





where Kp is a constant particular to a given type of explosive. Equations (3-5) and (3-6) are 
the classical "semi-empirical" expressions for the scaling relationships for the first oscillation 
period and maximum radius of an explosion gas bubble. 
B. NUMERICAL ANALYSIS MODELS 

Two different numerical analyses were conducted for deep explosion gas bubbles in 
the absence of a nearby boundary surface. The first analysis was designed to correspond to 
one of a series of experiments conducted by Swift and Decius (1950) at the Woods Hole 
Oceanographic Institution during and shortly after the Second World War. These researchers 
conducted a large number of experiments to determine constants for the semi-empirical 
equations (3-5) and (3-6) above for different types of explosives. The tests were conducted 
at a deep depth in deep water, to minimize the effects of buoyancy and boundary surfaces. 


The problem geometry for this analysis is shown in Figure 3-1. A very small (0.299 
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Figure 3-1. Geometry of Problem for Free-Field Bubble Analysis 


kg TNT equivalent) charge was detonated at a depth of 178.6 m in seawater. The 0.299 kg 
TNT equivalent charge weight was calculated by Swift and Decius to account for the extra 
energy of the detonator and booster material used above the energy of the TNT in the main 
charge. This particular shot (shot G72F) was chosen from all of the shots conducted because 
the bubble oscillation period and maximum bubble radius from this experiment were closest 
to the mean values from all of the experiments. Although the experiments used cylindrical 
charges with a height to diameter ratio slightly larger than one, a spherical charge was 
modeled for this analysis. 

In this experiment the maximum radius of the bubble was much less than the depth at 
which it was located. Thus, in the numerical analysis model for this problem, gravity was 
neglected and the pressure was assumed to be uniform in all directions from the charge. With 


this assumption, the problem is spherically symmetric. The seawater surrounding the charge 
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was treated as an invicid, irrotational fluid, and heat and mass transfer were assumed to be 
negligible over the time frame of the analysis. Herring (1950) has shown that the transfer of 
heat over one bubble oscillation cycle through conduction at the surface of the bubble 
constitutes a negligible fraction of the total energy of the explosion. 

Since this problem as modeled has spherical symmetry, a one-dimensional model is 
appropriate. In order to model this one-dimensional problem with the three-dimensional 
analysis code MSC/DYTRAN, a tall, thin pyramid shaped region of fluid was used. This 
accounts for the increasing volume as one moves away from the charge. Boundary conditions 
on the sides of this fluid volume were made "no flow" conditions, as flow across the 
boundaries is precluded by the flow in adjacent fluid volumes. The analysis model used for 
this problem is shown in Figure 3-2. 

At the apex of the pyramid shaped volume, a rectangular parallelpiped was used, with 
the equivalent volume a pyramid shaped volume would have. This avoids the use of a five 
sided pyramid shaped volume at the point, which is not a standard finite element shape. 
Rectangular parallelpiped shaped finite volumes were used to model the charge, and non- 
rectangular hexahedron finite volumes were used to model the remaining volume. Because 
the primary interest in this analysis was in the bubble behavior rather than the primary shock 
wave, only three finite volumes were used to model the charge. 

In order to prevent reflection from the far end of the finite volume of fluid modeled 
from affecting the results during the duration of the analysis, a very large volume of 
surrounding fluid was modeled. A total of 999 finite volumes were used for this analysis 


model, with the radial length of individual elements increasing with increasing distance from 
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Figure 3-2. Numerical Analysis Model for 0.299 kg TNT Explosive Charge 
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the charge. 

The TNT for this problem was modeled using a JWL state equation, with state 
equation parameters taken from the Lawrence Livermore National Laboratory Explosives 
Handbook (Dobratz, 1981). With this state equation, the pressure in the "burned fraction" 
of a charge is related to the specific internal energy and density by 


R, R, 

2 Oni, a Ls 
= -—— |e +B) 1-—— lex +@ E 3-7 
p | = | 4 p | p NP, (3-7) 


1 2 


where N=P/Pp, Po is the initial density, Eis the specific internal energy (per unit mass), and A, 
B, o, R,, and R, are constant parameters for the explosive. The "burned fraction" is just that 
portion of the explosive contained within a spherical detonation front traveling outward with 
detonation velocity d, which is another parameter characterizing the explosive. For TNT 
compressed to an initial density of py of 1.630 g/cm’, the JWL state equation parameters are 
(Dobratz, 1981): 

E=4.29 x 10° J/kg 

A =3.712 x 10" Pa 

B=3.231x 10’ Pa 

o = 0.30 

R, = 4.15 

R, = 0.95 

d= 6930 m/s 

In order to model the seawater in which this experiment was conducted, a polynomial 


state equation was used. This state equation relates the pressure in the fluid to its density p 


and specific internal energy £ per unit mass by 


2), 


p=a,p+a,n*+a,n°+(b,+b,n +b n”)p o£ (3-8) 


where 





(3-9) 


is the condensation and a@,, @, a3, 55, 5,, and 6, are constants for the fluid. Constant for this 
state equation for seawater, for use at condensation values pt of up to 0.8, were determined 
by fitting available Hugonoit data from the literature to this state equation form, and adjusting 
the bulk modulus and density to the accepted values for seawater. The resultant state 
equation parameters were (Chisum and Shin, 1995) 

a, = 2.306 x 10’ Pa 


a, = 8.432 x 10° Pa 
a, = 8.014 x 10’ Pa 


b, = 0.4934 
b, = 1.3937 
b, = 0.0000 


QO) = 1025 kg/m’ 

The initial conditions given to the seawater in this problem were an initial 
condensation value of zero (initial density of 1025 kg/m’), and an initial specific internal 
energy of 3750.4 J/kg. This initial specific internal energy was determined from equation (3- 
8); it represents the specific internal energy necessary to give the seawater an initial pressure 
equal to the hydrostatic plus atmospheric pressure for this problem. It was found necessary 
to use the initial specific internal energy rather than the condensation to set the initial 
pressure, because seawater is so incompressible that there is only a minuscule density change 


at this depth, and significant round off errors would be introduced by setting the initial 
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pressure with the condensation. 

The second model geometry analyzed was similar to Figure 3-1, except that an initial 
charge depth of 1000 m and an initial TNT charge of 3.5 kg was used. This corresponds to 
an initial charge radius of 8.0 cm. The objective of this analysis was to examine the flow 
characteristics in the surrounding fluid over one oscillation of the bubble. For this reason, 
uniform element spacing in the radial direction was used. A pyramid shaped volume of fluid 
was again used. In this case, the fluid was only modeled out to a distance of 16 m, and a 
simplified "non-reflecting" boundary condition was placed on the last element. This boundary 
condition was that the flow across the face of the element at the boundary would have the 
same values as in that element. 

The radial element spacing used was 0.125 cm, and a total of 12800 finite volumes 
were used. It was not necessary to use this many elements to get a "grid independent" 
solution, but the use of a fine mesh better illustrates the characteristics of the shock wave. 
The computer program used utilizes artificial bulk viscosity to control oscillations behind 
shock waves, and this smears the shock front over a number of elements. Mass, energy, and 
momentum are still conserved, but the shock wave doesn't look as steep. Thus, using more 
elements makes the shock wave "look" more like a shock wave, since it 1s smeared over a 
smaller distance. 

G NUMERICAL ANALYSIS RESULTS 

As stated in the last section, the first of the problems analyzed corresponded to an 

experiment conducted by Swift and Decius. Figure 3-3 shows the predicted radius time 


history from this analysis, and the first maximum bubble radius and period measured from the 
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Figure 3-3. Radius Time History for 0.299 kg TNT Charge Detonated at 178.6 m 
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experiment. This plot shows excellent agreement between the analysis results and the 
experimental measurement. Also shown on this plot, for comparison purposes, is the radius 
time history curve generated from step-by-step numerical integration of equation (3-3), which 
neglects both the internal energy of the gas bubble and fluid compressibility. This curve is 
seen to agree quite well with the numerical analysis curve up until the bubble reaches 
maximum volume, but to begin diverging from there. 

It is seen in Figure 3-3 that the numerical analysis predicts the bubble radius and 
oscillation period to decrease in subsequent oscillations. This is due to the radiation of a 
bubble pulse when the bubble collapses. This effect cannot be generated without including 
fluid compressibility. 

Figure 3-4 shows the pressure, impulse, and fluid particle velocity time histories, at 
a point located two maximum bubble radii from the center of the charge. This figure 
illustrates that, although the peak pressure in the first bubble pulse is much less than the peak 
shock wave pressure, the impulse generated by the bubble pulse is still large because of its 
longer duration. 

In their experiment, Swift and Decius also measured the second period and second 
maximum radius for this shot. Their measurements indicated that the second maximum radius 
was 11.6 inches (about 29.5 cm), and that the measured time for the second bubble minimum 
was 30.85 msec. These results indicate that an additional energy loss, above that due to 
acoustic radiation, occurred in this experiment, since acoustic radiation was accounted for in 
the numerical analysis for this problem. 


One possible explanation for this in that the bubble surface was unstable near the 
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Figure 3-4. Numerically Predicted Pressure, Impulse, and Velocity Time History for 
0.299 kg TNT Charge Detonated at 178.6 m, at a Point Two Maximum 
Bubble Radii From the Charge Center 
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minimum volume. Hicks (1970a) has noted that photographs of non-migrating bubbles 
typically show numerous needle-like water jets protruding into the bubble surface when the 
bubble is near minimum radius, and that this "spray" of water could significantly cool the hot 
bubble gases. 

In the second free-field numerical analysis, the flow external to the bubble was 
examined. Figure 3-5 shows the bubble radius time history for this analysis, and Figure 3-6 
shows the variation of the internal energy of the bubble as a function of time. The internal 
energy 1s seen to rapidly decrease from its initial value as the bubble expands, and then to 
increase moderately as the bubble is compressed by the surrounding fluid. 

Figures 3-7, 3-8, and 3-9 show the pressure in the fluid at 1 msec intervals. Because 
the shock wave is radiating spherically, the peak pressure of the shock wave decreases rapidly 
with increasing distance. The pressure behind the shock front drops below hydrostatic 
pressure, and this low pressure region extends for a considerable distance from the bubble by 
5 msec, when the bubble is near maximum radius. Also, by 5 msec there is a second sharp 
drop in the pressure behind the shock wave, and after 5 msec this sharp drop extends further 
and further into the fluid. 

This is caused by the reversal of flow in a portion of the fluid. The fluid immediately 
behind the shock front is still moving outward, but the fluid near the bubble surface is moving 
inward as the bubble contracts. As time goes on, fluid further and further away from the 
bubble begins moving inward. By the time the bubble reaches minimum volume, the pressure 
near the bubble is significantly higher than even the primary shock wave. As the bubble 


reaches minimum volume and begins expanding again, a moderately low pressure but long 
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Figure 3-5. Predicted Bubble Radius Time History for 3.5 kg TNT Charge Detonated at 
1000 m 
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Figure 3-6. Predicted Bubble Internal Energy Time History for 3.5 kg TNT Charge 
Detonated at 1000 m 
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Figure 3-7. Pressure Distribution From 1 msec Through 5 msec After Detonation of a 


3.5 kg TNT Charge at a Depth of 1000 m 
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Figure 3-8. 
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Figure 3-9. Pressure Distribution From 9 msec Through 13 msec After Detonation of a 


3.5 kg TNT Charge at a Depth of 1000 m 
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duration bubble pulse begins traveling outward away from the bubble. 

Table 3-1 compares the results between the experiment conducted by Swift and 
Decius and the current numerical analysis predictions. Extremely good agreement is seen for 
the first maximum bubble radius and period, while the second maximum bubble radius and 
period are appreciably less. The solution computed using the noncompressive theory 
neglecting the internal energy of the gas bubble was 8% lower than the actual period. Internal 


gas energy and fluid compressibility thus have an observable effect on the results. 





Table 3-1. Comparison Between Experimental and Numerical Analysis Results For First and 
Second Maximum Bubble Radius and First and Second Bubble Period. 
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IV. EXPLOSION GAS BUBBLE BEHAVIOR NEAR SIMPLE BOUNDARIES 


A. INTRODUCTION 

While free-field explosion gas bubbles dynamic behavior has been discussed in the 
previous chapter, practical problems typically involve the interaction of these bubbles with 
boundary surfaces. As discussed in the introduction, additional phenomena are observed 
when an explosion gas bubble is generated near a boundary surface. These include a change 
in the oscillation period of the bubble, and migration of the bubble from its initial position in 
the fluid medium. 

These phenomena will affect the damage caused by whipping of the hull girder of a 
naval vessel. The time interval between arrival of the primary shock wave and the first bubble 
pulse can be very near or equal to the fundamental period of the whipping mode of vibration 
of the hull girder, significantly increasing the whipping mode response of the vessel. Since 
the bubble pulse is generated at the conclusion of the contraction phase of the bubble 
oscillation, the period of this oscillation is obviously important. The migration of the bubble 
towards or away from a boundary during this oscillation is also important, as this determines 
the origination point of the bubble pulse. This affects both the arrival time of the bubble pulse 
and its magnitude, as the magnitude decreases with increasing distance. 

Since the available analytical models are known to not give quantitatively accurate 
results for the effects of boundary surfaces on underwater explosion gas bubbles, a number 
of numerical analyses using Eulerian finite volume meshes were conducted. These included 


a free-field analysis in which no nearby boundary was present, a series of analyses in which 
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Figure 4-1. Geometry of Problem for Bubble Near Boundary Analyses 


a plane infinite rigid boundary was located at vanous standoff distances from the charge, and 
a series of analyses in which a plane infinite free (constant pressure) surface was located at 
various standoff distances from the charge. 
B. NUMERICAL ANALYSIS MODELS 

The basic problem geometry for these analyses consists of a 10.24 kg cylinder of 
TNT, with a height to diameter ratio of one, located near an infinite plane boundary at a depth 
of 1000 m, as shown in Figure 4-1. Underwater explosion experiments also typically use 
cylindrical charges, for practical reasons. As in the free-field analyses described in the 
previous chapter, the choice of a deep charge depth serves to significantly reduce the period 
of oscillation of the bubble, and hence the analysis time. 

A deep charge depth also serves to simplify the analyses, as at a deep depth the 


variation in hydrostatic pressure in the surrounding fluid for a given depth change 1s 
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proportionately much smaller than at a shallow depth. The hydrostatic pressure in the fluid 
near a deep charge can thus be approximated as having a uniform value, the hydrostatic 
pressure at the depth of the charge, provided that the overall dimensions of the explosion gas 
bubble is much smaller than the depth of the charge. This was the assumption made for the 
analyses described in this chapter, for which the maximum bubble radius was less than 0.1% 
of the depth of the charge. 

The migration due to gravity can be separated from the migration caused by the 
presence of plane boundary surfaces, for cases in which the boundary surfaces have no 
component normal to the direction of the gravitational acceleration. Furthermore, deep 
explosion gas bubbles are known to experience little vertical migration due to gravity (Hicks, 
1970a). Thus, gravity can be neglected except for its effect on the hydrostatic pressure in the 
surrounding fluid. This permits separation of the effects of gravity from the effects of a 
boundary surface, facilitating investigation of the influence of the boundary surface on the 
dynamic behavior of the bubble. For the same purpose, Blake and Gibson (1981) have 
conducted experiments using spark generated vapor bubbles with the entire experimental 
apparatus in free fall. 

For these analyses, the cylindrical charge is assumed to have its axis normal to the 
plane boundary surfaces. With this orientation, the problems analyzed are axially symmetric, 
with the symmetry axis being the axis of the charge. To analyze these axially symmetric 
geometries using the three-dimensional code MSC/DYTRAN, wedge shaped volumes 
comprising two degrees of arc were used. 


Figure 4-2 illustrates the overall shape of the volume of fluid modeled, for the free- 
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Figure 4-2. Overall Geometry of Analysis Model for Free-Field Analysis 
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field analysis in which no boundary surface is present. The vertical line oriented in the y 
direction in this model is the axisymmetric symmetry axis of the problem. The boundary 
conditions imposed on the two semi-circular areas, lying in the planes oriented at +/- 1 degree 
about the axisymmetric symmetry axis from the x-y plane, are "no-flow" (rigid wall) boundary 
conditions, as flow across these boundaries is precluded by the flow in adjacent wedge 
segments (not modeled or shown). 

As discussed in the section in the previous chapter describing free-field analysis 
results, the dynamic behavior of underwater explosion gas bubbles is influenced by the 
conditions in the surrounding fluid adjacent to and at some distance from the bubble. Thus, 
a large volume of surrounding fluid was modeled in these analyses. This acts both to 
incorporate the necessary fluid and to avoid having to specify particular boundary conditions 
at the remote boundary of the finite volume of modeled fluid. This remote boundary was 
simply made a "no-flow" (rigid wall) boundary, with the distance to the boundary being made 
so great that shock wave reflection from it could not affect the analyses during the time frame 
of the analyses. For the free-field model volume shown in Figure 4-2, this remote boundary 
is the curved area connecting the semi-circles rotated +/- 1 degree about the y-axis from the 
x-y plane. These semi-circles are located at a distance of 400 m from the detonation point; 
as discussed below, an increased element size is used far away from the center of the charge, 
so that relatively few extra elements are needed to model fluid very remote from the center 
of the charge. 

The finite volume meshes used for the analyses described in this chapter were 


developed to meet several goals. It was important that the results obtained be independent 
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of the mesh used to obtain them, 1.e. that the solution was "grid independent." This goal 
could be met simply by using a large number of equal sized elements (finite volumes) for the 
analyses. However, if carried to extremes this practice is in direct conflict with other goals 
that a practical model must meet; that is, that the model must be able to run using a 
reasonable amount of computer resources (such as memory and storage space for results), 
and that the model must run to completion in a reasonable amount of time. It was also 
considered important that the results be independent of the shape of the elements used in the 
area near the charge, even if the bubble migrated during the analyses. A considerable amount 
of effort was expended in order to meet all of these goals for the analyses described in this 
chapter. 

During the development of meshes for these analysis, 1t was eventually determined 
that a mesh consisting of three separate regions could meet these goals. From a two- - 
dimensional perspective, the central region of these meshes consists of a rectangular region 
with square elements. This configuration allows elimination of element shape as a factor, 
since if this element geometry can still produce spherically shaped bubbles and the bubbles 
remain within this region during the analyses, then the shape of the element was not important 
in determining the results. The second region is a transition from this square element region 
to a spherically diverging type of region, which is necessary if the total number of elements 
is to be kept to a reasonable value. The third region is the spherically diverging type region, 
in which the dimensions of the element increase with increasing distance from the center of 
the charge, which, again, serves to limit the number of elements to a reasonable value. 


The use of this spherically diverging region is supported and suggested by the results 
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from the previous chapter, where it is seen that characteristics of the flow in the surrounding 
fluid are changing much slower in the remote fluid than in the fluid near the charge. As the 
purpose of the analyses described in this chapter is to determine the dynamic behavior of the 
bubble rather than the primary shock wave, it does no harm if the primary shock wave is 
spread over a larger region. And since a reasonable number of elements is still used in the 
circumferential direction (in the free-field model, 34 elements over 45 degrees, i.e. about 1.3 
degrees per element), this element geometry should still be valid for the analyses in which a 
boundary surface is present, and there is some circumferential change in the flow field. 

The dimensions of the rectangular (square meshed) region used was several times the 
maximum bubble radius of the charge. This was found to be necessary to eliminate mesh 
reflection affects on the bubble radius-time history for the free-field analysis. This meshing 
in sufficient to capture the important effects which occur within a short distance of the charge, 
as seen in the previous chapter. 

The mesh transition region provides a transition from the square meshed central region 
to the spherically diverging outer region, and reduces the number of elements in the 
circumferential direction per unit circumferential distance. This reduction is important in 
limiting the total number of elements needed for the analyses. Figure 4-3 shows the geometry 
of the central region and transition region for the free-field analysis model; this is just an 
expanded view of Figure 4-2, in the region in which the charge is located. 

To insure a grid independent solution, free-field analyses were conducted for which 
the square element dimensions in the central region were 2.5 cm and 3.333 cm (with slight 


differences in the gross dimensions of the central region (247.5 cm and 250 cm width)), and 
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Figure 4-3. Geometry of Analysis Model for Free-Field Analysis in Area Near Charge 
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the bubble radius versus time curves were found to lie on top of each other. It was concluded 
from these results that the 2.5 cm square dimension solution was a grid independent solution 
for the free-field case. It was then assumed that this meshing would be adequate for analyses 
involving a nearby boundary surface, since the meshing for these analyses was kept similar 
to the free-field meshing. 

Figure 4-4 shows the gross finite volume analysis mesh from a two-dimensional 
perspective, for the free-field analysis. Figure 4-5 shows the central rectangular region of this 
mesh, and the transition from this central region to the outer (spherically diverging) region. 
The total number of finite volumes (elements) used for this analysis is 47340. The central 
(square meshed) region uses 19602 elements, the transition region uses 338 elements, and the 
outer region uses the remaining 27200 elements. In this outer region, the mesh is designed 
such that the dimensions of individual elements remain nearly equal as the mesh diverges from 
the central region. 

The maximum bubble radius predicted by the free-field analysis was about 70.65 cm, 
and the first bubble period from this analysis was about 14.47 msec. These values were used 
to assist in quantification of the results for analyses which included a boundary surface. For 
analyses involving a boundary surface, lengths are given in terms of the non-dimensional 
standoff distance *, where h* is the standoff distance from the center of the charge to the 
nearest point on the boundary in units of maximum free-field gas bubble radii. Times are 
given in terms of the non-dimensional period 7*, where 7* is the time in units of bubble free- 
field first oscillation periods. 


The free-field axisymmetric model was used as the starting point for other analysis 
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Figure 4-4. ini 
gure 4-4, Overall Finite Volume Analysis Mesh for Free-Field Analysis 
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Figure 4-5. Finite Volume Analysis Mesh for Free-Field Analysis in Area Near Charge 
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models involving a rigid wall boundary and a constant pressure surface. Analyses were 
conducted at standoff distances of h* ranging from 1.062 to 4.034, for both types of 
boundaries. These h* values correspond to standoff distances from the center of the charge 
to the boundary ranging from 75 cm to 285 cm, respectively. 

Table 4-1 summarizes the characteristics of the models used in these analyses. The 
total number of finite volumes (cells) used in these models ranged from 28680 (for the h* 
equal to 1.062 analyses) to 42708 (for the h* equal to 4.034 analyses). The models for these 
analyses were similar to the free-field model described above, but with additional elements 
present below what was an additional symmetry plane for the free-field model (bisecting the 


axis of the charge). 
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Table 4-1. Summary of Characteristics of Finite Volume Analysis Models 


Figure 4-6 shows the geometry of a typical model used in these analyses from a three- 
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Figure 4-6. Overall Model Geometry for h*=1.380 Analyses 
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dimensional perspective. This figure is for the analyses in which h* has a value of 1.380. 
Figure 4-7 shows an expanded view of the geometry of this model, near the central region. 
In these figures, the lower of the triangular shaped areas parallel to the z-x plane is the ngid 
wall or constant pressure boundary. The upper (interior) triangular shaped area bisects the 
center of the charge, and was the symmetry plane for the free-field analysis. 

The overall finite volume mesh used for the model geometry shown in Figures 4-6 and 
4-7 is shown in Figure 4-8, from a two-dimensional perspective. Figure 4-9 shows an 
expanded two-dimensional view of this model in the region near the charge. As in the free- 
field model, this model has a square meshed central region, which for this model is a 
rectangular wedge 247.5 cm wide by 345 cm high. This is bordered on two sides by a 
transition region connecting this square meshed area with the more coarsely meshed outer 
region. 

The lower boundary parallel to the z-x plane in Figures 4-6 through 4-9 is made either 
a rigid boundary or a constant pressure boundary, by imposition of boundary conditions on 
the faces of the elements which lie on this boundary. For the ngid wall analysis, the 
appropriate boundary condition is a "no-flow" condition. For the constant pressure boundary 
analysis, the pressure on the element faces which constitute this boundary is made equal to 
the initial hydrostatic pressure in the surrounding fluid, for all time. The other flow 
parameters (velocity, density, and specific internal energy) at the boundary then take on the 
values in the finite volumes in to or out of which material is flowing. This constant pressure 
boundary is only an approximation of an actual air-water interface, as it does not account for 


changes in the shape of the boundary (which is assumed to remain planar). 
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Figure 4-7. Model Geometry for h*=1.380 Analyses in Area Near Charge 
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Figure 4-8. Overall Finite Volume Analysis Mesh for h*=1.380 Analyses 
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Figure 4-9. Finite Volume Analysis Mesh for h*=1.380 Analyses in Area Near Charge 
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The finite volume models for the other analyses conducted at vanous other standoff 
distances are very similar to the model shown in Figures 4-6 through 4-9. These models vary 
only in the amount of fluid added between the plane bisecting the axis of the charge and the 
boundary plane. 

The 10.24 kg cylindrical TNT charge used in the simulations discussed in this chapter 
is modeled as a 20 cm high, 20 cm diameter charge with an initial density of 1.630 g/cm’. 
Properties for this charge were modeled using a JWL form state equation, with state equation 
parameters taken from the Lawrence Livermore National Laboratory Explosives Handbook 
(Dobratz, 1981). This is the same state equation described in Chapter III, and the same 
parameters given there for this state equation were used for the analyses described in this 
chapter. In Figures 4-8 and 4-9, all of the elements initially contain seawater, with the 
exception of the 32 elements initially containing the TNT charge. This charge is shown as 
shaded elements in Figure 4-10. This TNT is detonated using a spherical detonation wave 
traveling outward from the center of the charge (midway up the charge on the axisymmetric 
symmetry axis) at a constant velocity of 6930 m/s. 

The seawater in all of the remaining elements was modeled using the same polynomial 
State equation and state equation parameters described in Chapter III. Initial conditions in this 
seawater were an initial density of 1.025 g/cm’, and an initial pressure of 10.153 MPa 
(0.10153 Kbar), which is the total hydrostatic pressure (atmospheric plus seawater head). As 
discussed in Chapter III, this initial pressure was set by specifying a nonzero initial specific 
internal energy (20076 J/Kg) to avoid round off errors (even at a depth of 1000 m, there is 


only a minuscule density increase). Using a constant volume specific heat value of about 1 
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Figure 4-10. Locate of Charge in Finite Volume Analysis Mesh for h*=1.380 Analyses 
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cal/g C, this works out to a temperature change of only 7.8 C, for which the results should 
not be affected. 

Several simplifying assumptions were made in these numerical analyses. As in the 
simulation described in the last chapter, the seawater was assumed to behave as an invicid and 
irrotational fluid. Heat and mass transfer between the seawater and the explosion gas bubble 
were also assumed to be negligible over the time frame of the analyses. Herring (1950) has 
shown the peak stresses due to viscosity are negligible in comparison with hydrostatic 
pressure, even using Taylor's (1943) velocity field (which is known to significantly 
overestimate the bubble migration velocity), and that the transfer of heat over one bubble 
oscillation cycle through conduction at the surface of the bubble constitutes a negligible 
fraction of the total energy of the explosion. 

G: SIMULATION RESULTS 

1. The Effect of Boundaries on Bubble Period in a Compressible Fluid 

As discussed in the introduction, the period of the first bubble oscillation and the 
displacement of the bubble during this time interval are important factors in determining the 
amount of damage that the bubble can cause. Figure 4-11 shows the volume equivalent 
spherical radius (the radius of a spherical bubble with the same volume) time history, for 
analyses in which a rigid and a constant pressure boundary were located at a standoff distance 
of 1.380 maximum free-field radii. The use of the volume equivalent spherical radius as the 
ordinate scale in this figure is not meant to imply that bubbles were actually spherical when 
a boundary was present. This ordinate scale is simply (3 V/4z)'”, and was chosen merely for 


convenience in comparing results from the various analyses with the free-field analysis results. 
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Figure 4-11. Bubble Radius Time History for Free-Field and h*=1.380 Analyses 
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Time (msec) 


The free-field radius-time history is also shown in this graph for comparison purposes. The 
period of the bubble is obviously influenced appreciably by the presence of the nearby 
boundary, as 1s the maximum volume of the bubble. Similar bubble volume equivalent radius 
time history data was generated for the analyses conducted at the other initial standoff 
distances, and used to determine the individual bubble oscillation periods. 

By compiling the bubble period results from the numerical analyses conducted at 
different standoff distances, the effect of the standoff distance to each type of boundary can 
be assessed. This is done in Figure 4-12, in which the non-dimensional bubble period 7* 
(where 7* is the first bubble period in units of free-field bubble periods) is plotted against the 
inverse of the non-dimensional standoff distance h*. For convenience, the results from 
analyses with rigid wall boundaries are shown on the same graph with results from analyses 
with constant pressure boundaries. The inverse standoff distance with an ordinate value of 
0.0 represents the situation in which the boundary is located at infinity (1.e. the free-field 
case). Points to the nght of this represent non-dimensional inverse standoff distances to a 
rigid wall, and points to the left are non-dimensional inverse standoff distances to a constant 
pressure boundary. 

Also plotted in Figure 4-12, for comparison purposes, is the non-dimensional period 
versus the inverse of the non-dimensional standoff distance (7* vsersus (1/h*)) relationship 
predicted by an approximate analysis due to Herring (1950). In Herring's analysis, the fluid 
is assumed to be incompressible, the gas bubble is assumed to have negligible internal energy 
and remain nearly spherical, and the bubble is assumed to be fairly remote from the boundary 


(terms of order 1/h*’ and higher are neglected). Herring's formula for the bubble period can 
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Figure 4-12. Bubble Period Variation with Standoff Distance to Rigid and Constant 
Pressure Boundaries 
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be expressed as 


rr Ben) eal 
4h“ 

where X,,, 1s the average radius of the free-field bubble during the first oscillation, and R,,,. 
is the maximum free-field radius. In this equation, the upper (plus) sign is for a rigid wall 
boundary, and the lower (minus) sign is for a constant pressure boundary. The R,,, value is 
determined by integrating the free-field radius vs time curve over the first period, and then 
dividing by that period; here, the free-field radius vs time curve from the numerical analysis 
was used, and integrated numerically. 

It is interesting that the 7* versus 1/h* curve generated from the numerical analyses 
has about the same slope as that predicted by Herring's analysis, for the mgid boundary type. 
Herring's analysis neglected terms of order (1/h**) and higher, and assumed small bubble 
shape deviations from spherical, so it is apparent that equation (4-1) could be in error for 
charges fairly near a boundary surface. Since Herring neglected compressibility of the 
surrounding fluid, and changes in the surrounding flow field caused by passage of the initial 
shock wave (which, as seen in the results in the last chapter, can extend over an appreciable 
distance from the bubble), equation (4-1) can also be shown to be in error for charges very 
remote from a boundary. In a compressible fluid, there is a finite standoff distance beyond 
which the absence or presence of a boundary can have no affect upon the first oscillation 
period of the bubble, owing to the finite wave speed in a compressible medium. Thus, in 
Figure 4-12, for a compressible fluid there must exist a finite region of abscissa values on 


either side of the 0.0 1/h* value for which 7* must be equal to 1.0. This distance can be 
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estimated, assuming a constant wave speed of about 1500 m/s in seawater, as half the distance 
traveled by shock wave in 14.47 msec (the free-field bubble period), i.e. about 10.85 meters. 
This corresponds to an h* value of about 15.4 (and a I/h* value of about 0.065). 
Furthermore, if the boundary is located somewhat closer it should have still have little effect 
upon the bubble period, as most of the oscillation of the bubble will have already occurred 
before the reflection of the primary shock wave gets back to the bubble. 

Thus, although equation (4-1) is the classical correction for the effect of ngid 
boundaries and free-surfaces on the oscillation period of an explosion gas bubble, a better 
correction can be made by including fluid compressibility. One question that arises is how the 
results from this investigation this might be applied underwater explosion gas bubbles 
resulting from other charge types, weights, and depths. One solution is to use a simplified fit 
through the data points in Figure 4-12. One proposed fit is a "bilinear" one through the data 
points for either type of boundary. The term "bilinear" is meant to imply that, for either a 
ngid wall or a constant pressure boundary individually, the curve consists of two linear 
segments (one of which has zero slope). This proposed fit is shown as a solid line in Figure 
4-12, and labeled as a "Bilinear Approximation." One limitation of this fit is that its accuracy 
decreases rapidly for standoff distances of less than two maximum free-field radu, for the 
constant pressure type of boundary (i.e. for 1/h* values greater than 0.5). For the ngid wall 
type of boundary, it appears to be reasonably accurate over the entire range of standoff 
distances shown in Figure 4-12. 

The advantage of this bilinear approximation lies in the relative ease with which it can 


be utilized for other charge type, weight, and depth parameters, as will be discussed shortly. 
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The functional form for this proposed peniod correction factor for the effect of a boundary 


surface on the first period of an explosion gas bubble is given by the expression - 





] Rive Ra 
T*=l+u a (4-2) 
h, 4h, 

where uw is a unit step function (zero if its argument is negative, one if its argument is 


positive), FR, 1s the average free-field bubble radius over the first period, X,,,. is the maximum 


free-field bubble radius, and h,’ is defined by 








ae (4-3) 


where /,* is an effective maximum standoff distance (in units of maximum free-field radii) for 
which the period of the bubble starts to be effected appreciably by the presence of a boundary 
surface. In Figure 4-12, this is the point at which the slope of the bilinear approximation curve 
changes from 0.00 to 0.20. For the bilinear approximation curve shown in Figure 4-12, this 
effective maximum standoff distance 1s given by 1/h)* = 0.27, 1.e. 4)* has a value of about 3.7 
maximum free-field radii. 

In order to apply Hernng's correction factor, equation (4-1), it was necessary to know 
the quantity R,... This implies knowledge of the free-field radius time history curve, but since 
the shape of this curve vanes little with the particular explosion parameters, the value of #,,, 
will generally be close to 0.8 R,,. To apply equation (4-2), it is necessary to know, in 


addition, the value of the parameter /,*. In general, this parameter may be a function of the 
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charge type, weight, and depth. However, it should be relatively easy to determine an 
empirical expression for this parameter for different types of explosives, since for a particular 
charge weight and depth, this parameter can be determined by a single accurate measurement 
of how much Herring's correction factor is in error for a single (well chosen) standoff 
distance. 

An example should clarify this. Suppose that an experiment were conducted using a 
particular charge type, weight, and depth, at a standoff distance of h*=2.00 from a rigid 
boundary, and Herring's correction factor was determined to predict a period increase of 
10%. The value of R,,/R,,,2 15 then 0.8, from equation (4-1). Ifthe actual period increase for 
this experiment was measured to be 4.6%, then, from equation (4-2), 1/h;=0.23. Thus, from 
equation (4-3), 1/h,* = 0.27. 

Again, the range of validity of equation (4-2) should be kept in mind. From Figure 
4-12, an appropriate limitation would appear to be /*>2.0 for a constant pressure boundary, 
and A*>1.0 for a rigid boundary. 

For calculation of the first bubble period at standoff distances less than two maximum 
free field bubble radu from a constant pressure boundary, a better approximation is needed. 
A "Shifted Quadratic" fit for the constant pressure boundary data shown in Figure 4-12 is 


defined by 





rt : a (4-4) 


* * 
TCP (hy cp) 


where wz is a unit step function and 
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_ (4-5) 


This characterization depends upon two parameters. The fit shown did not require a linear 
term in the shifted non-dimensional inverse standoff distance. The values for the shifted 
quadratic fit shown in Figure 4-12 fitting the numerical analysis data were m,-p=0.4 and 
h’,cp=0.17. The range of validity for the fit defined by equations (4-4) and (4-5) is h*>1, i.e. 
it is valid for any standoff distance from a constant pressure boundary. The fact that the non- 
dimensional period change seems to vary linearly with the non-dimensional inverse standoff 
distance for a rigid wall boundary and as the square of the non-dimensional inverse standoff 
distance for a constant pressure boundary may well have some physical cause or explanation. 

Ds The Effect of Boundaries on Bubble Migration in a Compressible Fluid 

The bubble center of mass displacement time history for the analyses in which a ngid 
boundary and a constant pressure boundary were located at an initial non-dimensional 
standoff distance of h*=1.380 is shown in Figure 4-13. This is the same case for which the 
bubble volume equivalent radius time history curves were shown in Figure 4-11. In Figure 
4-13, the displacement is taken to be positive if it is towards the boundary, and negative if it 
is away from the boundary. For the rigid wall boundary analysis, this figure shows a small 
initial migration of the bubble away from the boundary at early times (when the bubble is 
expanding). At later times (when the bubble is contracting) there is a significant migration 
towards the boundary. This migration is much larger than the initial migration away from the 


boundary, and the rate of migration increases as the bubble contracts to its minimum volume. 
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Figure 4-13. Bubble Center of Mass Displacement Time History for h*=1.380 Analyses 


69 


The opposite affect is seen when the boundary is a constant pressure boundary; the 
bubble initially moves towards the boundary, and later moves rapidly away from the 
boundary. The influence of fluid compressibility is also seen in Figures 4-11 and 4-13. There 
is an initial interval of time during which neither the type nor even the presence of the 
boundary affects the volume of the bubble, and it is not migrating. 

The displacement time histories, whose overall characteristic is described above, also 
show a small amplitude higher frequency component. The period of this small amplitude 
component ts about twice the time interval between the start of the analyses and the time 
when the bubble begins migrating, which seems to be roughly proportional to the standoff 
distance. Hence, this higher frequency component would appear to result from the reflection 
of the primary shock wave from the boundary. 

The general shape of the bubble displacement curves, as seen in Figure 4-13, has been 
qualitatively explained by Cole (1948, pp. 331-332): 


In the case of a rigid surface, the presence of the boundary interferes with 
radial fluid flow of water, whether outward or inward, near a spherical surface in 
its vicinity. Initially, when the pressure in the gas is in excess of the hydrostatic 
pressure, the water on the side of the bubble surface near the wall is less readily 
displaced, and the bubble surface moves away from the wall. The effect is 
relatively small, however, because the net pressure (in excess of hydrostatic) is 
positive for a short part of the bubble period, and the bubble is small during this 
time. When the pressure falls below hydrostatic, acceleration of flow toward the 
bubble surface does not occur as readily on the side toward the wall, and the flow 
must be such as to bring the surface nearer to the wall. A considerable amount 
of momentum is imparted to a large mass of water in this way when the bubble 1s 
large. As the bubble contracts, the momentum acquired becomes concentrated 
in a smaller mass of water near the bubble, and the velocity of flow in this region 
increases. The bubble surface must then move toward the wall with increasing 
speed as if attracted to it. This effect is so much larger than the repulsion when 
the pressure exceeds hydrostatic that the dominant motion is an apparent 
attraction increasing the bubble velocity toward the wall as it contracts, even 
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though the momentum of the flow is decreasing in the most contracted stages. 
A free surface has the opposite effect on bubble migration, as in this case 

the water at the surface is free to move but must due so in such a way as to 

equalize the pressure with that of the atmosphere. 

Hernng (1950) conducted an approximate analysis for the migration of bubbles near plane 
rigid and free surfaces, in which the method of images is employed for motion of a bubble in an 
incompressible fluid, the internal energy of the bubble is neglected, the bubble is assumed to be 
not too near the boundary surface (terms of order 1/h** and higher are neglected) and remain 
nearly spherical, and migration of the bubble is treated as a small correction to the motion of the 
bubble. Herring's formula for the migration velocity of the center of the bubble is 

gdh 31? dr 3 | “) - (4-6) 
dt Aprdt 2p2r7~° dt 
where ¢ is the time, / is the standoff distance between the boundary and the center of the bubble, 
and r is the radius of the bubble. In this equation the plus sign is for a ngid boundary, and the 
minus sign is for a free surface. 

Together with a knowledge of the function r(t), equation (4-6) can be separated and 
integrated to give an equation for the displacement of the center of the bubble. However, results 
obtained from doing so are known to be in poor agreement with experimental results, significantly 
over predicting the migration of the bubble late in the collapse phase (Cole, 1948, p.348). This 
is when the assumptions used by Herring in his analysis are least accurate, as the migration of the 
bubble at this time is not a small perturbation of the free-field bubble motion, and the time rate 
of change of the bubble radius is over predicted if internal energy is neglected. 


Nevertheless, it is worth noting that the first term on the right hand side of equation (4-6) 
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is a periodic term which changes sign when dr/dt changes sign, i.e. when the bubble is at its 
maximum radius. The second term is monotonic, and builds up to a large value when the 
magnitude of dr/dt becomes large. This is in qualitative agreement with the behavior seen in the 
displacement versus time results for the various numerical analyses. 

As discussed earlier, the location of an explosion gas bubble at the end of its first 
oscillation is important because this is where the bubble pulse originates from. Taking a cue from 
Figure 4-12, values for the bubble displacement at the time the bubble reached minimum volume 
(which vanes depending upon the boundary type and initial standoff distance) were compiled for 
the different numerical analyses, and plotted as a function of the inverse standoff distance. This 
plot is shown in Figure 4-14, in which the abscissa values represent the non-dimensional inverse 
standoff distance, with values to the nght of 0.0 being for a rigid boundary and values to the left 
of this point being for a constant pressure surface, as in Figure 4-12. The total displacements at 
the time of the first minumum is plotted in units of maximum free-field radii, so that a value of 0.4, 
for example, represents a center of mass displacement of 0.4*70.65 cm from the initial center of 
the charge. 

The effect of a finite wave speed in a compressible fluid is again seen in Figure 4-14; 
remote boundaries cause little or no displacement of the bubble. This plot shows a general 
linearity for bubbles fairly near either type of boundary. The only exception is the case in which 
the bubble was initially closest to a constant pressure boundary, for which the initial migration 
towards the boundary and the increase in the bubble radius over the free-field value cause the 
bubble to actually contact the boundary when the bubble was near maximum radius. Neglecting 


this left most point in Figure 4-14, the results show that the first period displacement near either 
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Figure 4-14. Bubble C.M. Displacement at First Minimum Variation with Standoff Distance 
to Rigid and Constant Pressure Boundaries 


ig 


type of surface can be well approximated by a bilinear curve when plotted against the inverse 
standoff distance. This curve is shown as a dashed line in Figure 4-14. The slope of the non-zero 
portion of this curve is about the same for either type of boundary, again suggesting a genera! 
semi-empirical formula for bubble migration near either type of boundary. 


The proposed functional form for this formula is 











| My 
; (4-7) 


where AY* is the displacement at the first minimum in free-field radii, w is a unit step function, my, 
is the slope of the non-zero slope portion of the "bilinear approximation" curve, as shown in 


Figure 4-14, and h,* is a effective standoff distance for bubble migration, given by 


a (4-8) 


where /* is the actual standoff distance in maximum free-field bubble radu, and #, * characterizes 
a maximum effective standoff distance for which appreciable bubble migration occurs. 

For the bilinear approximation curve shown in Figure 4-14, m, and h,* have the values 
of about 0.87 and 6.2, respectively. Again, these parameters might vary with charge type, weight, 
and depth, but can be determined empirically without great difficulty. The proposed bilinear fit 
given by equations (4-7) and (4-8) appears, from Figure 4-14, to be a reasonably good 
approximation for total bubble migration at the time of the first minimum for standoff distance 
values h*>1.25 for the constant pressure type boundary, and h*>1 for the ngid wall type 


boundary. 
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3. Bubble Shape Departure From Spherical 

While the time at which the bubble reaches its minimum volume and the displacement of 
the bubble at this time, and the effect of standoff distance on these quantities, were the primary 
items of interest in these analyses, the overall shape variation of an explosion gas bubble 
generated near a boundary surface is also of interest. In approximate analytical derivations such 
as for equation (4-1), the bubble is assumed to remain nearly spherical. It is of some interest to 
see how nearly true this assumption 1s, if only to form some opinion as to the validity of the 
derivation. While the analysis code used in these numerical analyses does not actually keep track 
of the location of the boundary surface between different materials within the mesh, it is possible 
to get an overall view of the shape of the bubble by examining the density of the various mesh 
volumes. 

Figure 4-15 shows the results for the analysis case in which a plane constant pressure 
boundary was initially located at a standoff distance of 1.380 maximum free-field radii from the 
center of the charge. In this figure the dark areas represent the bubble, and the times are given 
in terms of the first period of the bubble. The boundary plane is that plane below which times are 
indicated. The bubble is seen to be fairly spherical when near its maximum volume, but to have 
a pronounced "kidney shape" by time 80% of its first oscillation period has elapsed. It then 
continues deforming, eventually assuming a "spherical cap" shape by the time it reaches minimum 
volume. This figure also shows the migration of the center of mass of the bubble away from the 
boundary. 

Figure 4-16 shows the variation of the bubble shape with time when a rigid boundary was 


located at 1.380 maximum free-field radii from the center of the charge. Again, the dark area 
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Figure 4-15. Bubble Shape Vanation with Time for Bubble Located at h*=1.380 from a 
Constant Pressure Boundary 
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Figure 4-16. Bubble Shape Variation with Time for Bubble Located at h*=1.380 from a 
Rigid Boundary 


represents the bubble, and the boundary plane is that plane below which times are indicated. In 
this figure, bubble shapes are only shown for the last 5% of the first oscillation period; at 95% of 
the first bubble oscillation period the bubble is still nearly spherical. By 98% of the first 
oscillation period, the bubble has become "kidney"-shaped. 

From these figures, it might be expected that derivations assuming bubble sphericity 
would be better at this standoff distance for the ngid boundary case than for the constant pressure 
boundary case. Examining the bilinear approximation in Figure 4-12, which can be thought of 
as a compressibility corrected approximation to Herring's derivation of equation (4-1), this would 


seem to be the case. 
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V. CONCLUSION 


The current understanding of the effects of fluid compressibility on the dynamic behavior 
of underwater explosion gas bubbles is very limited; the existing theories for underwater 
explosion gas bubble phenomena neglect entirely the effects of fluid compressibility. AS a 
consequence, these theories fail to accurately predict important phenomena associated underwater 
explosion gas bubbles. It was for this reason that the research described in this dissertation was 
undertaken. 

In the first part of this investigation, it was shown that a finite volume based numerical 
analysis technique which incorporated fluid compressibility could accurately account for the 
behavior of an explosion gas bubble in the free-field case up until the time of the first bubble 
minimum, and it was verified that the radiation of energy in the bubble pulses does not account 
for all of the energy loss that occurs between bubble oscillation cycles. 

In the second part of this investigation, the behavior of underwater explosion gas bubbles 
near rigid and constant pressure boundary surfaces was investigated, including effects arising 
because of fluid compressibility. The failure to account for fluid compressibility in the current 
theoretical development has made predictions for bubble behavior based upon this theory more 
qualitatively correct than quantitatively accurate. It was found possible, when fluid 
compressibility is taken into account, to quantitatively describe the behavior of explosion gas 
bubbles quite close to a boundary surface. 

By conducting analyses at a number of different standoff distances from these boundaries, 


it was discovered that first bubble oscillation period and the migration of the bubble at the time 


i 


it reached minimum volume could be quantitatively described over all standoff distances of 
interest when fluid compressibility was taken into account. Functional relationships between 
standoff distance and the first bubble period and the net bubble migration were developed which 
have a form very suitable for semi-empirically scaling the results of this investigation to situations 
involving other charge types, weights, and depths. 

It is recommended that a series of carefully designed experiments be conducted to both 
verify the formulas developed in this investigation and to develop the neccessary empirical 
constants used in these formulas, so that the effects of boundary surfaces can be better taken into 
account in the future. These results of this investigation should form a foundation for improving 
the accuracy of calculations for the whipping response of ships and shipboard equipment due to 


detonation of nearby mines or torpedoes. 
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